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ABSTRACT 


The  paper  is  the  third  and  final  part  in  the  series  of  three  devoted 
to  the  detailed  analysis  of  the  three  basic  versions  of  the  finite  element 
method  in  one  dimension.  The  first  part  [1]  analyzed  the  p- version,^ the 
second  part  [2]  concentrated  on  the  h  and  h-p  version* and  the  present 
third  part  addresses  the  adaptive  h-p  version. 


Accession  For  / 


1 .  INTRODUCTION 

This  paper  is  the  third  and  final  part  in  the  series  of  three  which 
analyzes  the  h,  p  and  h-p  versions  of  the  finite  element  method  in  one 
dimensional  setting.  It  has  been  shown  in  Part  1  and  2  that  the  selection 
of  the  mesh  and  degree  of  elements  is  essential  for  the  performance  of  the 
method.  We  have  shown  that  the  proper  selection  of  the  h-p  version  leads 
to  the  exponential  rate  of  convergence  while  the  h-version  with  improper 
mesh,  e.g.  uniform  mesh,  gives  very  low  algebraic  rate  when  a  singularity 
is  present.  The  adaptive  approaches  are  essential  for  solving  complex 
problems,  because  the  structure  of  the  solution  is  not  known  a-priori. 

In  recent  years  the  adaptive  methods  came  to  be  in  the  focus  of 
interest.  Various  papers,  see  e.g.  [3]  [4]  [5],  address  the  question  of 
adaptive  approaches  in  the  Finite  Element  Method.  In  two  dimensional 
adaptive  research  code  FEARS  (see  [6])  and  PLTMG  (see  [7])  are  avail¬ 
able.  Both  codes  deal  with  the  h-version  and  linear  (p  *  1)  elements. 
However,  there  is  no  adaptive  h-p  version  code  and  only  little  work  has 
been  done  addressing  this  question.  See  {8]  [9], 

In  this  paper  we  analyzers  theoretical  frame  of  the  adaptive  h-p 

version  and  based  on  it  we  provide  concrete  algorithm  for  the  one  dimen- 

-  •.! 

sional  problem.  It  is  proven  that  in  the  case  that  the  solution  has  x®- 
type  singularity,  the  adaptive  algorithm  give  an  exponential  rate  of  con¬ 
vergence,  very  close  to  the  optimal  one  analyzed  in  the  second  part  of  the 
paper. 

We  expect  that  the  principles  used  here  in  the  one  dimensional  set¬ 
ting  could  be  successfully  applied  also  in  the  higher  dimensional  case. 
Although  above  we  used  the  notion  of  an  adaptive  approach  in  broad  sense, 
we  will  later  distinguish  in  a  more  precise  setting  between  feedback  and 


adaptive  approach  (cf.  [10]  [11]  [12]  [13]).  By  feedback  approach  we 
understand  the  approach  when  previous  (computed)  information  are  sequen¬ 
tially  used.  The  adaptive  approach  is  in  a  feedback  which  has  well 
defined  optimal  properties.  The  distinction  between  the  feedback  and 
adaptive  approach  is  often  worthwhile  in  a  more  precise  analysis. 

We  will  develop  in  this  paper  an  abstract  frame  of  the  adaptive 
approach  and  its  theory.  We  are  concerned  here  only  with  the  convergence 
and  its  rate  in  the  energy  norm.  Section  2  focusses  on  the  algorithm. 
Section  3  deals  with  its  convergence,  and  Section  4  analyzes  its  rate  of 
convergence.  Section  5  gives  some  numerical  and  a  short  discussion  of 
implementational  aspects.  Section  6  summarizes  the  major  properties  of 
the  three  basic  versions  of  the  FEM. 

2.  THE  ABSTRACT  SETTING  OF  THE  ADAPTIVE  h-p  VERSON  ALGORITHM  OF  THE 

FINITE  ELEMENT  METHOD. 

First,  we  will  make  some  definitions.  A  mesh  A  is  a  partition  of 
interval  [0,1].  For  convenience,  a  mesh  A  may  be  regarded  as  a  set  of 
nodal  points  or  a  set  of  non-overlapping  closed  intervals,  the  union  of 
which  is  [0,1].  The  number  of  intervals  contained  in  a  mesh  A  is 
called  the  cardinality  of  a  mesh  is  denoted  by  m(A).  To  each  mesh  inter¬ 
val  1^  we  assign  a  positive  integer  p^,  which  is  called  the  degree  of 
the  mesh  interval.  These  degrees  constitute  the  degree  vector  The 

superscript  A  indicates  its  relation  to  the  mesh  A;  if  there  is  no 
confusion  it  will  be  often  omitted. 

Definition  1.  S  is  the  set  of  mesh-degree  combinations.  Its  element 


is  called  the  pair,  where  A  is  a  mesh  and  j>  is  its  associated  degree 
vector. 

E  is  also  regarded  as  a  set  of  the  pairs  (I,p)  where  I  €  A 
and  p  is  the  degree  of  I. 

Definition  2.  Let  E  €  S,  the  number  of  degrees  of  freedom  of  E  _is_ 

m(  A) 

(2.1)  N  =  deg(E)  =  l  pj. 

i-1 


We  can  make  S  a  partially  ordered  set  by  defining  the  following 
partial  ordering  on  5. 

1  2 

Definition  3.  Let  E^,  E^  €  S.  E^  =  (A ,£  )  ,  ^2  ~  t^ien 

E^  <  E^  if  and  only  if: 

1)  A2  is  a  refinement  of  A^  i.e.,  as  sets  of  nodal  points  one 
has 


Ai  S  &2 


and  we  will  write 


41  '  42 ’ 


2)  £  is  a  refinement  of  £  ,  i.e.,  A^ 


<  A^ ,  and  if 


▼  (  1  )  ,  A 

li  6  V 


I(2)  €  A 
1  "  A2  ’ 


l(2)  c  -j-O) 

Li  "  j 


then  p£L)  <  p^2). 


this  case  we  write 


-CD 


(2) 


,  (1)  ,  (2) 
(p.  and  p' 


In 


are  the  degrees  associated  to  the  intervals 


We  will  call  £2  a  refinement  of  £ 


Definition  4.  The  local  error  function 


E( [  a  ,b] ,p) 


is  a  non-negative  real  valued  function  defined  on  the  set 


{0<a<b<l,  p  >  1,  pis  integer}. 


which  satisfies  the  following  hypotheses: 

(El)  E([a,b],p)  is  continuous  in  a,  b.  It  is  non-increasing 
in  a  and  p,  and  non-decreasing  in  b. 

(E2)  E([a,a],p)  =  0  for  0  <  a  <  1,  p  >  1. 

(E3)  (p-approximability) .  For  any  fixed  [a,b]  £  [0,1] 


lim  E( [a ,b] ,p)  *  0. 

p-HJO 

For  some  given  X,  1  <  X  <  <»: 

(E4)  (reverse  sub-additivity).  If  eg  [a,b],  p  >  1,  then 


(E([a,c],p)X  +  E([c,b],p)X}1/X  <  E([a,b],p). 


OO 

(E5)  (h-approximability) .  Let  {A^}^  be  a  sequence  of  meshes 


for  which 


1 im  max  {  1 1 j }  =  0 

n+»  I€A 


lim  {  £  E(1, 1  )X)  ^X 

n-*-®  Ie  A 


0. 


(if  x  =  «, 


a.  >  0,  we  define  {  £  a  }  =  max  a.). 

1  i=l  l<i<m  1 

The  number  X  is  called  the  index  of  the  local  error 
function. 

Remark  1 .  These  hypotheses  are  very  natural  if  we  consider  the  error 
function  to  be  the  local  error  of  the  best  approximation  in  a  certain 
function  space.  For  example,  let  u  £  L^COjl)  and 

E([a,b],p)  =  inf  llu-vll  , 

v£P  ,  L2(a,b) 

p-1 

then  it  is  easy  to  check  that  all  hypotheses  are  satisfied  for  X  =  2. 

We  have  the  following  simple  corollaries  which  follow  immediately 
from  the  definition: 

Corollary  1.  If  Ij  £  I2 ,  then 


(2.2) 

Ed^p)  <  E(I2,p) 

(VP  >  !)• 

Corollary  2. 

If  a  = 

x0<x1<***<xk=b,  p  >  1,  then 

(2.3) 

k 

{  l 

j  —  1 

E(  txi_i  jXf]  ,p)*>  <  E([a,b],p). 

i=l 


Proof.  (2.2)  follows  from  (El).  (2.3)  follows  from  (E4)  by  induction.  I 


Definition  5.  A  local  error  indicator  e([a,b],p)  is  a  non-negative  real 
valued  function  defined  on  the  same  set  as  E([a,b],p).  And  there  is  a 


constant  0  <  C  <  1,  independent  of 


Ce( [a,b]  ,p)  <  E([a,bj,p)  <  C-1e( [a ,b] ,p) 


(2.4) 
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It  is  obvious  that  for  the  local  error  indicate  we  have: 


Corollary  3.  The  following  properties  hold: 


(E2)' 

e( [a, a] ,p)  =  0  for 

all  0  < 

a  <  1 ,  p  >  1 

(E3r 

lim  e(  [a,b]  ,p)  =*  0 

p+oo 

for  any 

[a,b]  £  [0,1  ] 

(E5)' 

Let  (A  }  ,  be  a 

n  n*l 

sequence 

of  meshes,  1 

and 


lim  max  {  1 1 ( }  =  0 

n+®  I£A 

n 

then 

lim  {  l  e(I,l)V/X  -  0. 

n*00  I  £  A 

n 


(X  is  the  index  of  local  error  function). 


Definition  6.  The  global  error  based  on  the  pair  E  =  (A ,j>)  is  given  by 

(2.5)  E  CD  =  {  l  E(I.,p.)V/X 

Ij€  A 

where  1  <  X  <  «,  px  is  the  degree  of  1^. 

Similarly,  the  global  error  estimator  based  on  the  pair  E  is  given 
by 

(2.6)  e.(E)  =  {  l  eU.,p.)V/X. 

A  I.€A  1  1 

i 

Clearly  we  have 

Corollary  4.  Let  C  be  given  (2.4),  1  <  X  <  «,  E  be  any  pair,  then 

(2.7)  Cex(E)  <  EX(E)  <  C_1ex(E). 

We  now  define  the  feedback  h-p  version  algorithm  (we  will  call  it 


the  algorithm  below) 


First,  let  the  local  error  indicator  e([a,b],p)  be  given.  We  will 
divide  the  intervals  of  any  mesh  into  two  categories,  called  the  tv- 
intervals  and  the  p-intervals .  We  will  also  say  that  an  interval  is 
of  h-type  or  p-type.  The  type  of  an  interval  is  defined  by 

Definition  7.  Let  0  <  y  <  »  be  given.  Let  E  =  (A,j>)  be  a  given 
pair,  I  6  A  and  p  the  degree  of  I,  let 

R  *  (R  -  0  if  e(I,p)  =  0). 

e( I ,p) 


Then  if  R  >  y,  I  is  said  to  be  an  h-interval;  if  R  <  y,  I  is  said  to 

be  a  p-interval  (with  respect  to  y  and  E). 

The  number  y  is  called  the  type-parameter.  Usually,  we  are 

interested  in  the  case  0  <  y  <  1. 

The  feedback  algorithm  is  now  defined  in  a  recurrent  way: 

Let  0  <  9  <  1  be  a  given  number,  called  the  refinement-parameter 

E  =  (A  ,2°),  and 
n  n  -c- 


»(n) 

max 


,,(n)  (n). 

max  e(I.  ,p.  ) 

l<i<m(A  ) 
n 


with  (I^n\p^n^)  €  E  .  (The  interval  (  1  on  which  e^n^  is 

l  i  n  j  n  max 

realized  will  be  called  the  critical  interval.)  Denoting 


eU) 


.d'n>.p'n)) 


then  all  intervals  I 


(n) 


€  A  with  e^n^  <  0 
n  i 


(n) 


e"  '  will  stay  the  same 
max 


(n) 


in  the  mesh  An+^  of  the  new  pair  En+^  and  the  degrees  p£  '  will  be 


also  kept.  If  for  some  interval  €  A_  on  which  e.^n^  >  9  •  e^n\ 


9 


then  there  are  two  cases: 


1)  if  is  an  h-interval,  then  it  will  be  bisected,  and  the 


degree  ^  will  be  inherited  by  both  the  bisected  intervals; 

2)  If  1^  is  a  p-interval,  then  remains  an  interval  of 


n+1 


but  its  degree  is  assigned  to  be  p^n^  +  1 


We  write 


-  T<V  s  T(In,e,-„6) 

where  e  stands  for  the  local  error  indicator,  y  the  type-parameter 
and  0  the  refinement-parameter. 

Definition  8.  The  above  rule  T  of  constructing  the  pair  from  an 

existing  pair  is  called  the  transition  operator  of  the  algorithm. 

The  subset  {E^} n=0  of  where  En  a  T C E ^  >  for  n  =  1,2,...,  is 

called  a  traj ectory  of  the  transition  operator. 

It  is  obvious  that  we  have 

Corollary  5.  A  trajectory  {E^}  is  a  monotone  increasing  sequence  in 
3,  namely,  E^  <  En+1  for  all  n  =  0,1,2 . 

Remark  2.  There  are  two  degenerated  cases  for  the  algorithm:  if  y  = 

0,  then  all  intervals  are  of  h-type  and  this  algorithm  gives  a  feedback 
h-version.  If  y  =  <*>,  then  all  intervals  are  of  p-types  and  will  never 
be  bisected;  in  this  case  we  will  obtain  a  feedback  p-version. 

Remark  3.  We  are  speaking  about  feedback  algorithm  because  the  current 
information  steers  the  flow  of  algorithm.  Often  such  algorithm  is  called 
also  adaptive  (see  e.g.  [10],  pp,  49-50,  [11]).  We  shall  distinguish 
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between  a  feedback  algorithm  and  an  adaptive  algorithm  in  the  sense  that 
the  adaptive  algorithm  is  a  feedback  having  well  defined  optimality 
properties  (see  [12],  [13],  [4]).  In  the  next  section  we  will  prove  that 
the  algorithm  is  convergent  and  hence  it  is  adaptive  with  respect  to  the 
convergence  measure. 

3.  THE  CONVERGENCE  OF  THE  ALGORITHM 

Definition  9.  If  for  any  trajectory  {^n}™=0  °f  the  transition  operator 

lim  E  (E  }  =  0, 

X  n 

n+<*> 

then  the  algorithm  is  said  to  be  convergent. 

We  will  prove  the  feedback  algorithm  defined  above  is  convergent. 
First  we  observe  that  the  (E5)  implies: 

00 

Lemma  1 .  Let  {iTn} n_Q  be  a  sequence  of  sets  of  non-overlapping  closed 
intervals  (not  necessarily  covering  the  entire  interval  [0,1]).  If 

lim  max  {  1 1 1 }  =  0, 

n+«°  I€r 

n 

then 

lim  {  l  E(I,1)V/X  =  0. 

n+"  L€  7t 

n 

Proof.  There  is  A  such  that 

-  n 

1)  as  sets  on  non-overlapping  closed  intervals, 

TT  c  A  . 
n  -  n 

2)  max  {  1 1 1 }  =  max  { 1 1 ]  } . 

I€A  I€w 

n  n 


we  obtain  using  (E5) 


Applying  (E5)  to  the  mesh  sequence 


<o 


n;  na0  * 


{  l  E(I,1)X}1/X  < 

I€ir 

n 


{  l  E(I,1)X}1/X  *  0 

l€An 


as  n  -*•  ®. 


As  consequences,  it  is  easy  to  obtain: 

(  n) , ® 

Corollary  6.  If  {I  ^n-Q  a  secluence  of  closed  intervals, 

p^n^  >  1  and  lim  j I ^ |  =  0,  then 
n+® 

lira  e(I(n),p(n))  =  0. 

n->® 


Corollary  7, 


Let  (tt  }  be  given  as  in  Lemma  1. 
n  n=u  ° 


(n) 


Assign  1^  6  7rn  an  integer 


(n) 


p^  >1  and  assume  that 


then 


1 im  max 


n-*-®  T(n) 
li 


lim  {  l 


n+®  ,(n) 

1  .  €  TT 

i  n 


/T(n)  (n) .A. 1/A 

I ^  >  P^  )  J 


=  0, 


We  now  prove 


Theorem  1.  The  feedback  algorithm  is  convergent. 

00 

Proof.  Let  {I  }  .  be  a  trajectory  of  the  transition  operator  J. 

-  n  n=0 


Corollary  5 


First  we  show  that 


(3.2) 


lim  max  {E(I^n^ ,p^n^ ) }  *  0 

n+«  x(n)€A 


where  An  is  the  mesh  of  £n,  and  p^n^  is  the  degree  of  (  a  , 

By  (2.4)  it  is  enough  to  show  that 


(3.3) 


i*  f  /  T  ( ( n)  v . 

lim  max  {e(I  ,p  )} 
n>«  r(n)€A 


Suppose  this  is  not  true,  then  there  is  a  subsequence  {nk}  and  a 


number  e  >  0  such  that 


(3.4) 


(O 

I  €A 


(n.)  (n  ) 

{e(I  k  ,p  k  )>  > 


(n.)  (nk) 

Furthermore,  let  I„  *  a  be  the  critical  intervals,  i.e.  e  is 

(nk)  K 

realized  on  1^  .  Corollary  6  implies  that  there  is  another  subsequence 

of  {n,},  we  still  denote  it  by  {nk} ,  such  that 


(nk) 

for  some  h  >  0,  for  otherwise  e  ■+•  0.  Since  I  <  E  ,  two 
.  raax  nk  nk+l 

^nk  'V+l 

intervals  I  and  I  are  either  non-overlapping  or 

c  c 

(n,  )  (nt+i)  (n,  ) 

I  <  z>  x  K  1  .  Because  x  k  c  [0,1],  there  can  only  be  a  finite 
c  “  c  c 

number  of  non-overlapping  intervals  with  length  >  h  >  0.  Thus  we 


conclude  there  is  again  a  subsequence  of  {nt},  we  again  denote  it  by 


{n,J,  such  chat  for  some  k  >  0 


T(nk+1)  D  (nk+2^ 

I  -  I  J  •  •  • 

c  ~  c 


^nk+l^  ^nk^ 

Each  I  is  either  a  result  of  several  bisections  of  I  »  or 

c  c 


T(Vi>  _  /V 

C  C  ( 
k'  such  that  I 


Because  by  our  assumption  |I 


>  h,  there  is 


is  never  bisected  for  k  >  k'  and  hence  I. 


the  p-interval.  Because  Ic  are  critical  intervals,  we  conclude  that 


(nk) 

lim  p  =  +  oo 

.  c 

k+°° 


(n  )  (n,  )  (i\) 

where  p  is  the  degree  of  I  .By  (E3)',  lim  emax 

k>» 

(n,  )  (n,  ) 

lim  e(I  ,  P  )  «  0.  This  contradicts  (3.4). 

i  c  c 

k+°» 

For  X  *  «  theorem  follows  immediately. 

Assume  now  that  1  <  X  <  °°.  Since  (3.1)  implies 


W  »  W  >  W  > 


(El)  and  (E4),  it  suffices  to  show  that  for  each  e  >  0,  there  is 


N(e)  >  0  such  that 


EX(ZN(e))  < 


For  each  h  >  0  we  define 


Sn  -  Cl  €  U  0  A  ;  | I !  >  h} 

k*0  n=k  n 
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where  &n  are  regarded  as  sets  of  intervals.  Then  for  each  h  there 
is  Nfl(h)  >  0  such  that  if  n  >  N«(h),  then 


tt  =  A\Z 
n  n' 


contains  no  interval  I  with  ( I j  >  h.  Let  hj  +  0  and  choose 

corresponding  n^  =  NqOij)  to  be  such  that  nj  +  ®.  Then  the  sequence 
00 

{ tt n  }  .  q  satisfies  the  condition  of  Corollary  7  and  thus 

j  J" 


lim  I  l 
j+®  j  (n  ) 

(I.  J  ^ 

\  i  n. 

J 


(n.)  (n.)  X 

e(Ii  2  ,pi  J  ^ 


=  0. 


Therefore  there  is  Jq  *  j(e)  such  that 


(3.5) 


l 

(n.  ) 


(n.  )  (n.  )  , 

E(I,  °,p,  ) 


<  z,2in. 


Jo 

h.  h. 

If  Z  =  0,  the  proof  is  finished.  Suppose  Z  ^  0.  Let  M  be  the 

hj 

J0 

number  of  intervals  contained  in  Z  .By  the  first  part  of  the  proof, 


there  is  >  jg  such  that 


(n.  )  (n.  ) 

J 1  3 1  ] 

!(Ii  ,p.  )  <  e/(2M) 


)  h^  h. 

for  all  I  1  €  Z  0  (recall  that  z  0  contains  the  intervals  which 


are  never  bisected  againt).  Therefore  we  have 
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(3.6) 


(n  )  (n  ) 

J  1  J  1  * 

[  Ed,  ,P,  ) 

(n  )  h  1 

It  Si  0 


<  e/ 21/x. 


Furthermore,  the  intervals  contained  in  A  \  I  are  either  inherited 

n. 

J1 

from  ir  or  the  results  of  several  times  bisections  of  the  intervals  of 

n1 

30 

it  ,  (3.5)  and  (El),  (E4)  imply  that 

% 


(3.7)  (n  )  L 

h  l€V'zh. 

J1  J1 


and  (3.6),  (3.7)  give 


(n  )  (n  )  X 
J  j  J  i 

E(I.  ,P1  )  i 


<  c/21^. 


E.  CE  )  <  e. 

\  n . 


Thus  the  proof  is  complete. 


Remark.  4.  Since  in  the  proof  we  did  not  use  any  information  of  the 
parameter  y,  we  have  shown  a  convergence  of  all  feedback  h,  p  and 
h-p  version  algorithms. 


Remark  5.  We  introduced  a  family  of  feedback  algorithms  which  create 
trajectories  {2^}.  We  define  the  performance  measure  Uq  of  the 
algorithm  so  that  if  E^d^)  ♦  0  for  any  trajectory  {E^}  of  the  algo¬ 
rithm,  then  Ug  =  1,  otherwise  Uq  =  0;  and  we  define  the  optimal 
performance  measure  to  be  with  maximal  value  (1  in  this  case).  Thus  we 
can  say  that  our  feedback  algorithm  is  adaptive  with  respect  to  the  per¬ 
formance  measure  up  (called  convergence  measure,  see  [12],  pp.  7-8). 


LV-V.'-.V 


6 


4.  THE  RATE  OF  CONVERGENCE  OF  THE  ALGORITHM 

In  order  to  study  the  rate  of  convergence  of  the  algorithm  it  is 
necessary  to  have  more  knowledge  on  the  local  error  function.  Motivated 
by  the  results  of  Part  1  [1],  we  will  study  the  algorithm  on  the  class  of 
local  error  functions  which  satisfy  the  following  hypotheses: 

(Al)  There  is  a  point  C  €  (0,1],  called  the  singular  point.  If 
€  €  I  ^  [0,1]  ,  then 


e( i » l )  <  c0 f 1 1 CT, 

where  o  >  0,  Cq  >  0  are  constants  independent  of  I. 

(A2)  There  is  a  non-increasing  continuous  function  <j>  :  (0,®)  ♦  (0,1) 
with 

lira  <fr(t)  *  1 

t*0 

lira  *(t)  *  0 

t+ao 


such  that  for  any  e  >  0  ,  p  >  1,  I  £  [0,1],  5  4  I,  and 
t  »  *  c>  cheire  is  C(e)  >  0  such  that 

E( I , p)  <  C(e)[*(t)]p. 


(A3)  If  £  4  I,  I  -  (0,1], 

E( I , p)  >0,  p  >  1 

<r 

the  constant  0  <  <  <  1  is  independent  of  I,  and  p. 

Observe  that  the  hypotheses  (Al),  (A2)  are  the  characterisics  of  the 
best  L2“approximation  error  of  the  analytic  function  with  an  xa-tvpe 


then  there  is  r  =  r(I)  >  0,  such  that 


..  E(I,P+1) 

E(  I , p) 
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singularity  at  £ .  In  Part  1  [1]  we  have  shown  that  if  the  local  error 
function  E(I,p)  represents  the  I^-error,  then  a  =  a  +  y,  and 
4>(t)  =  1  / ( 1+2 t+2/ 1(  1  +  t)  ).  The  results  of  Part  1  also  show  that  for  the 
function  (x-£)^  the  hypothesis  (A3)  is  satisfied. 

Lemma  2.  Let  the  local  error  function  E(I,p)  satisfy  (Al)  ~  (A3), 
e(I,p)  be  the  local  error  indicator  of  E(I,p),  then  e(I,p)  also 
satisfies  (Al)  ~  (A3)  with  different  constants.  More  precisely,  we  have 

(Al)'  If  £  £  I  c  [0,1],  then 

e(I,l)  <  C^l I) a. 


where  a  is  the  same  as  in  (Al),  Cq  >  0  independent  of  I. 

( A2 ) '  If  £  (  I,  I  c  [0,1],  t  =  *  c  >  °*  then  there  is 

C'(e)  >  0  for  which 

e(I,p)  <  C'(e)  [4>(  t)  ]  p 

where  p  >  1,  <j>  is  the  same  as  in  (A2). 

(A3)'  If  £  4  I,  I  c  [0,1],  then  there  is  r'  =  r'(I)  >  0,  such  that 


<  '  r  ' 


(Ll2±12 

(I.p) 


r' 


with  p  >  1,  0  <  <'  <  1 ,  <'  is  independent  of  I  and  p. 


Proof . 
have 


(Al)',  (A2)'  are  obvious.  Suppose  (A3)  holds,  then  by  (2.4)  we 


e ( I , p+ 1 )  C  lE(l,p+l)  -2 

e(I,p)  CE(I,p) 


e(  I ,  p-t-1 )  CE(I,p+l)  f2 

ed.P)  C-‘E(I,p 
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where  r  =  r(I)  >  0,  0  <  C  <  1,  0  <  k  <  1.  Let 

r'  =  r'(I)  =  C_2r(l), 

,  4 

K  =  C  K, 

and  (A3)'  is  satisfied.  | 

By  the  hypothesis  (El)  it  is  easy  to  see  that  (Al)  thus  (Al)'  may  be 
extended  to  p  >  1 : 

Corollary  8.  If  £  €  I  £  [0,1],  p  >  1,  then 


E(I.p)  <  CQ| I | 


e(I,p)  <  C'|l| 


hold  uniformly  with  respect  to  1  and  p. 


Lemma  3 .  Let  E(I,p)  satisfy  (A2)  and  (A3).  Then  we  have 


(4.1) 


r  <  <  ^q 


where  r,  <  are  defined  in  (A3),  q  =  f(t)  is  given  in  (A2)  and  t  = 
dist(S.I) 

- re¬ 


proof.  Bv  (A3)  there  is  r  such  that 


«  ‘  »  ‘  '• 


E(I,p)  >  C(  <r) 


This  implies 
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with  p  >  1,  C  being  a  constant  independent  of  p.  Comparing  with  (A2) 
we  then  obtain  for  any  p  >  1 


with  t  =  — fI).  Thus 


<r  <  q  =  4>(  t) 


„  -1 
r  <.  k  q, 


Corollary  9.  If  E(l,p)  satisfies  (A2),  (A3),  then 


(4.2) 


r '  <  k'  ^q 


with  r"  ,  <'  given  in  Lemma  2,  q  =  <j>(t)  as  before. 

Proof .  Obvious.  § 

Lemma  4.  Let  E(I,p)  satisfy  (A2),  (A3),  0  <  y  <  1  is  the  type- 
parameter  (as  defined  in  Section  2),  and  6  >  0  be  determined  such  that 


<{>(<5)  =  <'y 


1 


where  0  <  y^  <  Y  <  1  and  <'  is  the  constant  in  (A3)'. 
10,1], 
x  ^  ,  and 


If  x  €  I  L 


(4.3) 

then 


I  < 


L*ziL 

1+5  ’ 


„  _  e(I,p+I) 

R  =  — ; ;  •  <  y 

e(I,p) 


where  p  >  1.  Therefore,  I  is  a  p-interval. 


Proof .  Recall  0  <  k'  <  l  (see  Lemma  2).  Since  0  <  y^  <  y  <  1,  there 
exists  5  >  0  satisfying  f>(5)  =  <'Yj, 


(4.3)  implies 
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dist(C.I)  .  |x-C|-  I 

m  '  iii  >  5‘ 


therefore  <J>(t)  <  <p(S)  =  k'y^,  and  by  (A3)'  and  (4.2) 


e( I ; p+1 )  /K-\  ^(K'y  ) 

e(I,p)  }  K  V 


Y1  <  Y  > 


thus  I  is  a  p-interval. 


We  now  study  how  the  algorithm  does  if  the  hypotheses  (Al)  -  (A3) 

)QQ 

^_q  starts 

with  the  mesh  Aq  =  {[0,1]},  and  =  (Pq)»  Po  >  1*  In  this  case,  any 
meshes  of  the  trajectory  can  only  contain  the  interval  of  the  form 


(4.4) 


[■—,  \]  ,  1  <  k  <  2n,  n  -  0,1,2,, 

2  2 


Definition  10.  An  interval  of  the  form  (4.4)  is  called  a  binary  interval 


of  level  n. 


Lemma  5 .  Let  the  local  error  function  satisfy  the  hypotheses  (A2)  -  (A3). 
Let  y  be  the  type-parameter  0  <  y  <  1  •  If  x  (  [0,1],  x  ^  5,  then 
the  algorithm  will  generate  a  binary  interval  I  such  that  x  6  I,  I  is 
a  p-interval  and  it  remains  to  be  a  p-interval  in  the  further  process. 

Proof .  By  the  proof  of  Theorem  1  we  have  shown  (without  (A2)  -  (A3)) 


lim  max  {e(I^V^  ,p^ )} 

V+oo  (v)  ^  ^ 

I .  't  A 
J  v 


Therefore,  if  x  €  1^  €  Ay,  x  t  5  and  e(  I  ^ v)  ,Oj V ) )  >  0,  in  the 

further  process  either  will  be  bisected  or  will  go  up.  By 


Lemma  4  the  bisection  must  stop  if  the  size  of  the  interval  is  small 


enough,  in  other  words,  there  is  v.  such  that  if  x  (■  I 

0  '  j 


0' 


then 


I(V 

j 


is  a  p-interval  and  lies  in  all  A  .  v  >  . 

v  0 


Lemma  6 .  Let  the  local  error  function  satisfy  (A2)  -  (A3),  and  6  >  0  be 
defined  in  Lemma  4,  Suppose  that  the  pair  I  is  generated  by  the  algo¬ 
rithm  and  it  has  its  smallest  interval  which  is  of  level  n,  then  the 
total  number  of  intervals  in  the  mesh  A  of  £  is  bounded  by 


(4.5) 


M  =  2(  n+1 )  (L+l )  +  1 


where  L  =  [ 1+6 ] . 


Proof .  By  Lemma  4,  if  x  €  I ,  x  s*  £  and 


H  <  ^  |x-C|, 


then  1  is  permanently  a  p-interval  (i.e.,  it  remains  as  a  p-interval 
in  further  process).  This  is  true  if 


(4.6) 


H 


1 

1+L 


x-el 


where  L  =  [1+6], 

If  this  interval  is  of  level  k,  then  |l 


1 


and 


(4.7) 


Without  loss  of  generality  we  can  consider  the  interval  which  is 
right  to  Denoting  d  =  dist(£,I),  (4.7)  implies  that  for  each  k 


there  are  no  h-intervals  of  level  k  with 
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(4.8) 


d  > 


1+L 


Because  each  p-interval,  except  [0,1],  is  obtained  by  bisecting 
an  h-interval,  the  above  implies  that  there  are  no  p-intervals  of 
level  k  with 


(4.9) 


d  > 


1+L 

2  k- 1  £  k 


2L+3 
„k  ‘ 


mos 


Now  within  a  distance  from  £  between 
•1+L 


1+L 

r.k+1 


and 


1+L 
„k  ’ 


one  can  at 


t  put  in  [— —  h-intervals  (of  level  k)  or  (1+L)  p-intervals. 


Either  way,  the  number  of  intervals  within  above  range  is  bounded  by 

(1+L).  This  is  true  for  all  0  <  k  <  n.  Within  a  distance  ranging  from 
1+L 


0  to 


there  can  be  at  most  (1+L)  intervals  since  the  smallest 


interval  is  of  level  n.  Hence  the  total  number  of  intervals  to  one  side 
of  Z  will  not  exceed  (n+l)(l+L).  Including  the  interval  containing 
C ,  total  number  of  intervals  in  the  mesh  then  will  not  exceed 

M  =  2( 1+L) (n+1)  +1.  ® 


Lemma  7 .  Let  the  local  error  function  satisfy  (Al)  -  (A3),  0  <  y  <  1 , 

a o 

0  <  9  <  1,  and  let  {Z  }  q  be  a  trajectory  starting  with 


v  = 
"0 


([0,1] ,P0),  Pq  >  1.  In  addition,  let  the  local  error  indicator 
e(I,p)  satisfy  the  hypothesis: 


(El)'  if  I.lplj  -  [0,1],  IL  .  I2,  p  >  1,  then 


e(Ij.p)  <  e(I2,p) 


e(  I , p+1 )  <  e(  I , p) 


Then  we  have 

deg( E  ) 

(4.10)  E  (Z  )  <  C,[m(A  )]i/Ae  V 

where  Ay  is  the  mesh  of  Iy,  m(Av)  is  the  number  of  intervals  in 
A^,  deg(Z^)  is  given  by  (2.1),  and  Cj,  C£  are  positive  constants 
independent  of  v. 

Proof .  First  we  claim  that  there  is  a  constant  C  independent  of  v  such 
that  if  (l,p)  €  Zy,  then 

(4.11)  e(I,p)  <  Cyp. 

We  prove  this  by  induction  on  p  >  pg.  Note  that  the  trajectory  is  an 
increasing  sequence  of  pairs 

(4.12)  r0  <  <  r2  <  —  . 

Let  I j  p j  denote  an  interval  which  has  degree  p.  Suppose  for  some  v, 
X[p]  €  Av’  and  p  *  p0*  Let 

Mpoj  -  IIp0+U  -  •••  -  Itp-1 1  -  ^p] 

be  the  sequence  of  the  successsive  ancestor  intervals  of  Ij  j ,  each  of 
which  has  corresponding  degree  as  indicated  by  the  subscript.  By  the 
hypothesis  (El)'  we  have 

<*•13)  e(Ir  ,,p  )  <  e([0,l],p  )  =  CYP°. 

0  0  0 

We  not  let  p  >  pg.  Suppose 
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I 


e(I(p_1],p-0  <  CyP*1. 


If  ^P) 


I [ p_ i ] »  then  1 [ p_ i ]  is  a  p-interval. 


Thus 


e(l[p],p)  <  Ye(I[p],p-l) 

=  YeCljp.j ] ,P~1)  *  CTP. 

If  I[p]  I[p_! ] ,  then  there  is  Ifp]  -  Ifp_n  such  that  I(p]  £ 
xtp]  C-  I[p-l]  thuS  by  (E1)' 

e( 1 1  p] ’ ?)  <  e<I*[p]»P) 

<  Ye(l *[p_1],p-l)  <  Cyp. 


This  proves  (4.11). 

Let  e^  be  the  maximal  local  error  indicator  for 
max 

« ( v ) 

0  <  9  <  1  we  define  p  by 


E  .  For 
v 


(4.14) 


9e 


(v-1) 

max 


v  —  1,2,3,... 


where  C  is  given  in  (4.11). 


By  the  hypothesis  (El)'  and  (4.12)  it  is 


easy  to  see  that 

e 

(v?) 

>  e  , 

max 

max 

and  it  follows  immediately 

that 

/V 

(4.15) 

P 

<  P 

We  now  claim  that 


(4.16) 


P  <  p  , 
*max  K  ’ 


for  all  v  >  1 


( v) 

where  p  is  the  maximal  degree  of  the  intervals  of  E  .  We  will  prove 
max  v  K 

this  by  induction.  Observe  that  (4.11)  and  (4.14)  imply  that 


CYP°  >  «(10,11,P0)  > 


;<»-i 


=  cY 


n  „  -CD  , 

Pq  <  P  -  1. 


Since  p  equals  either  Pq  or  pg  +  1  we  obtain 


(1)  „  «(1) 

P  <  p 
rmax 


Suppose  we  have 


(v-1)  »(v-l) 

Pm-,„  <  P 


Now  for  p  ,  either  p 

L  m  o  v  '  L  n 


(v-1) 


or  p 

1  77 


(v-1) 


+  1.  In  the  first 


(v)  «(v-l)  *(v) 

<  p  <  p  ; 


and  in  the  second  case  the  interval  I  having  degree  p  was  a  p- 

max 

interval  of  E  , .  Hence  according  to  the  algorithm  we  have 
v-1 

,T  (v-1)  .  (v-1)  p^-i 

e(I,p  )  >  0e  =  Cy 

max  max 


By  (4.11)  this  gives 


Therefore  (4.16)  is  true  for  all  v  >  1.  From  (4.14)  and  (4,16)  we  obtain 


Then 


>) 

max 


< 


(v-1) 

e 

max 


P(V) 

rmax 

y 


deg  ZN 
C  m(  A) 

-  •  Y  • 

0Y 


m(  A  ) 
v 


A  AX 
v  n 


I  e(l.\Pp 

i*l 


1/X 


< 


C  ,  / »  vvl/X 

W  m  Av  } 


deg  l 

v 

m(  A  ) 
v 
Y 


and  (4.10)  follows. 


In  particular 


We  now  prove  the  main  theorem. 

Theorem  2.  Let  hypotheses  (Al)  -  (A3)  and  (El)  hold,  0  <  Y  ^  1> 

0  <  0  <  1.  Let  {£v)”=q  fee  a  trajectory  starting  with  =  ([0,1], p0 ) > 
p  >  1.  Then  there  exist  positive  constants  C1  and  C2  independent 
of  v  such  that 

Vo 

-C0(a  deg  £y)  - 

(4.17)  Exav)  <  Cf 


f where  1  <  X  <  00 


is  the  index  of  error 


function  (see  Definition  4)  and 


a  is  the  exponent  in  hypothesis  (Al)). 


Proof .  First  if  m(Ay)  is  bounded  by  a  finite  number,  then  by  Lemma  7  we 


will  have  an  estimate  C^e 


-C2deg 


which  is  better  than  (4.17). 


Therefore  we  can  assume  m(A  )  +<*>  for  v  -*■  ®. 

v 


Suppose  that  the  smallest  interval  I  of  X  has  a  level  n, 


thus  III  =  2  n.  Let  be  the  maximal  local  error  indicator. 

max 


Furthermore,  let  J  be  the  parent  interval  from  which  I  was  obtained  by 


bisecting  J  in  Ey  ,  <  v  -  1.  Therefore  J  is  an  h-interval  of 


level  n  -  1,  and  by  (4.8)  we  must  have 


(4.18) 


dist(£,J)  < 


L+l 


,n-l 


where  L  is  given  in  Lemma  6. 

a  A  A 

Let  I  be  such  an  interval  that  5  €  I  and  J  £  I,  it  is  easy  to 

A 

obtain  I  with 


ti|  <  ^±r+  1 


L+2 


2n”l  2n_1 


,n-l 


By  (Al)',  we  obtain 


Ki, p)  <  c0| 1 1 CT  <  c0(pjrf)CT 


where  p  is  the  degree  of  J  in  X^  ,  and  by  (El)'  we  get 


e(  J ,  p)  <  e(!,p)  <  * 


Because 


J  was  bisected,  we  must  have 


Using  (El)'  we  derive 


(4.19) 


•>>  < 

max 


<V  c 

e  < 
max 


9  ^e(J,p) 


where  C 


9_1C0(2L+4)°. 


On  the  other  hand  Lemma  6  shows  that 


m(A  )  <  2(n+l)(L+l)  +  1 

v 


because  the  smallest  interval  of  A  is  of  level  n.  Thus 

v 


m(A  )-l 
n  *  20+0“ 


and  (4.19)  shows  that 


ra( Ay)-1 

,<»>  <  c(— 

max  2° 


-  I  .  v-  m(A  ) 

-  c'<±>  a  >  v 


for  some  constant  C'  >  0,  independent  of  v.  Therefore 


(4.20) 


a  In  2  (  .  , 

/ y  ,  .  ..l/\  2(L+l)  m(V 

e  (Z  )  <z  C  (m(A  ))  e 
X  v  v 


In  Lemma  7  it  has  been  shown  that 


(4.21) 


e.  (Z  ) 
X  v 


< 


C(m(Av) 


1/X 


-(In 

e 


T  m(  A  ) 


(4.20)  and  (4.21)  give 
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e  (E  )  <  C(m(A  ))  e 

A  V  v 


lra  In  2  ,  .  , .  „  1  deg(Zv\ 

1/x  “  2 ^ 2(L+1 )  m(Av)  ln  Y  m(Av)  J 


<  C(m(Av))1/A  e 


/in  2  In  1/y  / - ; — / -  - 

-  / "Ta+T) —  '  /0  deg(V 


Noting  that  m(A^)  <  deg(A^) ,  the  above  inequality  implied  (4.17)  by 

taking  ^2  ~  ^  ^(L-t-l)'  ^  “  e  (£  small  enough  so  that  >  0),  and 

then  choosing  Cj.  £ 

In  (Al)  we  assume  £  €  [0,1].  Suppose  £  f  [0,1]  but  (A2),  (A3) 
hold,  then  we  have 

Theorem  3.  Let  hypotheses  (A2)  -  (A3),  (El)'  hold  with  £  |  [0,1],  0  < 

Y  <1,  0  <  0  <  1,  and  let  be  a  trajectory  starting  with 

=  ([0,1] ,pQ) ,  p0  >  1,  then  the  number  of  intervals  m(Ay)  is  bounded 
by  a  finite  number  when  v  -*  «.  Therefore 


-C  deg(E  ) 

(4.22)  EX(Zv)  <  C!e 

with  Cj ,  C2  >  0  independent  of  v. 

Proof .  By  Lemma  5,  since  in  this  case  x  /  £,  the  algorithm  will  generate 
a  permanent  p-interval  containing  each  x  £  [0,1]  (these  intervals  will 
never  be  bisected  again).  Clearly  there  are  only  finitely  many  such  inter¬ 
vals.  The  rest  of  the  part  of  the  theorem  follows  from  Lemma  7.  | 


We  now  discuss  the  adaptivity  of  the  algorithm.  Recall  in  Part  2  we 


obtained  for  the  model  problem  that  the  optimal  rate  of  convergence  in  the 
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energy  norm  of  arbitrary  mesh-degree  combinations  (the  pairs)  when  £ 
(0,lj  is  bounded  below  by 

C(a)  - - -  [(/I  -  1 ) 2  ]  ^  a"  l/2  )N 

/jf  0_  V2 

where  N  is  the  number  of  degree 
local  error  function  was  given  by 

E([a,b],p)  = 

and  A  =  2,  a  =  a  -  V2  .  We  see  that  this  rate  of  convergence  is  of  the 

-C  /oN  _  _2 

form  C^e  “  with  C2  =  £n[(/2  -  1)  ]. 

It  can  be  shown  under  the  assumptions  of  Theorem  2  (with  certain 

condition  on  <p )  that  there  are  constants  ,  C2  >  0 

-C2/oN 

(4.23)  $(a,N)  =C1e 

is  the  best  possible  estimate.  Therefore  we  define  a  performance  measure 

u  as  follows: 
e 

If  there  are  constants  C  >  0,  p  >  1,  such  that 
E.  (Z.  )  <  C[*(a,deg(Z  ))]l/p 

A  J  V 

for  v  =  0,1,2,...,  then  ue  =  1,  otherwise  ue  =  0.  We  then  can  say: 


of  freedom.  This  is  the  case  that  the 


inf  He'll.  ,  ,  s 

,  L2U'b) 
P~1 


Theorem  4.  Under  the  conditions  of  Theorem  2,  the  algorithm  is  adaptive 
with  respect  to  the  performance  measure  u  . 
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Remark  6.  The  notion  of  the  adaptivity  (as  the  optimality  of  the  feedback) 
can  be  defined  in  various  ways  here.  It  directly  relates  to  the  question 
of  comparison  of  feedback  algorithms.  In  [10]  [11]  it  is  shown  that  if  one 
considers  only  the  worst  case  problem  from  a  class  F,  then  for  many 
classes  F  the  nonadaptive  approach  is  as  good  as  any  adaptive  one. 

In  [12]  [13]  [4]  the  optimality  is  defined  asymptotically  (for  high 
accuracy  and  the  performance  of  a  trajectory  created  by  an  algorithm  (for 
every  particular  problem)  is  compared  with  the  performance  of  the  best 
trajectory  (for  a  given  particular  problem).  In  [4]  the  set  F  of  pro¬ 
blems  is  characterized  for  which  as  class  of  feedback  algorithm  create 
trajectories  with  comparable  performance  as  the  theoretically  best  trajec¬ 
tory  and  hence  the  feedback  algorithm  is  an  adaptive  one. 

In  this  paper  we  judge  the  algorithm  how  it  performs  with  respect  to 
a  worst  case  in  a  narrow  class  of  solutions  having  a  singularity  of  the 
type  xa  inside  or  outside  of  the  interval  I.  It  is  clear  that  not  a 
single  non-feedback  algorithm  can  perform  better  than  our  feedback  algo¬ 
rithm  for  this  class  of  solutions. 


Remark  7 .  It  is  possible  to  obtain  the  results  of  this  section  only  on  the 
base  of  the  hypotheses  (Al)'  -  (A3)'  and  (El)'  of  the  local  error  indicator 
without  the  assumptions  of  local  error  function  E(I,p).  In  fact,  the  only 
statement  in  the  section  which  required  the  property  of  E(l,p)  was 


(4.24) 


where  A 

v 


1  •  rT(v)  (v) 

lim  max  e(I  ,p.  ) 

i<v,ea 

i  v 


is  the  mesh  in  the  trajectory  {l^} 


=  0 


v=0  * 


Now  we  prove  this  directly  based  on  (Al)'  -  (A3)'  and  (El)'. 
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Let  be  Che  critical  interval,  and  be  c^e  middle 

point  of  I^V^*  Since  0  <  *  1*  there  is  a  subsequence  { v " }  5  {v} 

(v') 

lim  x  =  x  . 


Also  let  h  =  I  ,  since  0  <  h  <1,  we  can  assume  for  the  same  sub- 

V  1  C  V 


sequence  {v'} 


lim  =  h. 

V  ^  -*<» 


If  h  >  0,  there  can  be  only  finitely  many  1^  which  are  different. 

(  \)'  )  A 

Therefore,  there  is  v-,  if  v  >  cn,  then  I  =  I  will  be  fixed  as 
’  0  O’  c  c 

(v') 

a  permanent  p-interval,  therefore  ?c  00  and  the  algorithm  gives 


-  .  <V  Pc(',')-p'U0) 

e(Ic,P;  ')  <  e(Ic,pc  )  •  Y  -  o 


as  v'  +  <*>.  By  (El)'  the  maximal  local  error  indicator  is  non-increasing, 
therefore  (4.29)  holds. 

We  now  assume  h  =  0.  First,  ohserve  in  this  case  we  cannot  have  xc 
4  ?.  Form  (A2)'  and  (A3)',  the  argument  in  proving  Lemma  4  shows  that  if 
an  interval  does  not  contain  %,  then  it  cannot  be  bisected  infinitely 
many  times.  Thus  there  are  no  invarvals  with  dist(£,I^V^)  >  t  >  0 

and  !  1^  J  ^  j  -*■  0. 

This  shows  we  can  assume  xc  =  For  each  a  ,  let  the  inter¬ 
val  I^J  contain  both  %  and  ^ .  Clearly,  we  can  choose  ^  to 


be  such  that 
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I  *(v' ) 


15  -  x 


(v') 


+  h  , 
v 


( as 


0). 


Therefore  bv  (Al)'  and  (El) 


,T(v  )  (v  ).  /i(v  )  , 

e(I  ,p  )  <  e(I  ,1) 

c  c 


«  c0|!<“')!”  <■  o. 


This  completes  the  proof  of  (4.29). 


As  a  last  remark  we  indicate  that  the  above  results  are  all  valid  if 
we  increase  the  degree  uniformly  on  all  intervals.  In  fact,  the  proofs  are 
concerned  with  the  worst  possible  degree  distribution  made  by  our  algo¬ 
rithm.  If,  instead  only  increasing  the  degree  on  some  p-intervals  as 
described  in  Section  2,  we  increase  all  degrees  by  1  whenever  there  is 
p-interval  on  which  degree  is  supposed  to  be  increased  (cf.  Definition  8), 
then  we  will  obtain  an  adaptive  h-p  version  algorithm  which  produces  uni¬ 
form  degree  vector,  and  this  modification  does  not  cause  any  change  in  our 
original  proof. 

The  case  having  uniform  degree  vector  is  important  because  it  is  much 
easier  to  make  implementation  in  2  and  3-dimensional  case. 


5.  NUMERICAL  RESULTS 

Table  1  and  Figure  1  are  the  numerical  results  obtained  by  using  the 
adaptive  h-p  version  algorithm  described  in  Section  2.  The  problem  is  the 
model  problem 

-u"  =  f 

u ( 0 )  =  u(l)  =  0 
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with  the  solution 


u( x)  =  (x-5)+  -  (l-£)a*x  -  (-£)“(l-x) 


and  we  compute  for  the  case  a  =  0.7 ,  a  =  1.1,  5=0,  5=  0.3, 

respectively.  The  local  error  indicators  are  exact  local  error  of  the 
finite  element  solution.  According  to  the  theory,  we  will  obtain  an 
exponential  rate  of  error  reduction. 


(5.1) 


where  N  is  the  number  of  degrees  of  freedom  and  r  =  10. 

We  use  linear  regression  to  find  the  constants  C  and  x.  Comparing 
with  the  theoretic  values  <'  for  the  h-p  extension  with  geometric  mesh 
and  linear  degree  vector  when  qQ  =  0.1715  (the  optimal  one),  q  =  0.5 
with  corresponding  optimal  s,  we  obtain  the  following  table: 


TABLE  1 


a 

c 

< 

x//a-  V2 

<' 

<opt=1715 

q=0 . 5 

q-0.1715 

q=0 . 5 

0.7 

0 

2.458 

0.3174 

0.7097 

0.3424 

0.3036 

0.7656 

0.6789 

0.3 

6.280 

0.3566 

0.7974 

1.1 

0 

0.6106 

0.5563 

0.7181 

0.5930 

0.5259 

0.3 

6.115 

0.6642 

0.8575 

Figure  5.1  shows  that  the  error  reduction  curves  are  near  to  straight  lines 
in  the  /N  ~  loglell  scaled  graph,  as  expected.  The  slopes  shown  in  Fig. 
5.1  are  the  theoretical  ones  for  q0pt  and  q  =  0.5. 

For  the  implementation  of  the  algorithm,  as  mentioned  before,  the 
assumptions  of  the  theory  are  satisfied  for  the  model  problem,  -u"  =  f, 


when  the  solution  can  be  obtained  on  each  mesh  interval  separately.  In  the 
general  case,  the  local  errors  are  affected  by  global  error  and  the  assump¬ 
tions  are  not  to  be  satisfied.  Furthermore,  there  is  a  problem  in  finding 
an  effective  local  error  estimator  for  large  mesh  intervals  and  high 
degrees.  Although  there  are  difficulties  in  both  theoretical  and  practical 


IK)  50  100  200 

N 

Figure  5.1. 


aspects  of  the  adaptive  h-p  version  of  the  FEM,  a  program  was  written  for 
the  one-dimensional  two-point  boundary  value  problem: 

-(a(x)u')'  +  b(x)u  =  f ( x )  x  6  (0.1) 

u(0)  =  u(l)  =  0 

Figures  5.2(a),  (b)  are  the  results  obtained  by  using  the  adaptive 


h-p  version  FEM  program  to  solve  the  following  problem 
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with  the  solution 

u(x)  =  (x-£)+  “  (1-C)ax 


for  a  =  1.1,  E,  =  0  and  E,  =  0.3.  The  graph  is  in  /N  -  log II ell,,  scale, 

h, 

(a)  shows  the  error.  (b)  (the  dotted  curves)  are  curves  of  estimated  error 
(by  the  global  error  estimates). 

We  also  use  linear  regression  to  obtain  the  constant  k  defined  in 
formula  (5.1),  which  is  shown  on  the  figures. 

Our  results  show  the  program  also  performs  very  well  on  this 
problem.  The  global  error  estimator  is  very  reliable.  When  E,  =  0.3  the 
singular  point  will  never  be  a  nodal  point  during  bisection.  In  this  case 
the  curve,  although  oscillating,  still  gives  the  expected  exponential  rate, 
with  the  rate  of  convergence  and  the  error  itself  better  than  when  £  ^  0. 

The  program  basically  agrees  with  the  algorithm  as  described  in  Sec¬ 
tion  2.  There  are  nevertheless  some  different  features  in  the  program: 

1)  We  not  only  increase  the  polynomial  degrees  but  also  allow  lower¬ 
ing  them.  The  reason  for  doing  this  is  to  make  the  local  error  more  equi¬ 
librated  so  that  we  can  avoid  unnecessarily  increasing  the  total  number  of 
degrees  of  freedom.  This  is  done  in  the  following  way. 

For  each  v  =  1,2,...,  we  let  «  y  be  given, 


(5.2) 

with 


■  /  _(v). 

=  min{p  ,R  } 


v  =  1,2,, 


,(v) 


=  rain<RlR  =  ee(i^)}*  <I»P) 


V- 


Tf  /T(v)  (vK  , 

If  e .  =  e(  I .  ,p .  )  <  p 

3  3  3  v 


.  (v) 

degree 


g(v)  (e^V^  =  max  {e!V^}),  then  the 
max  maX  Uj<m(A  )  J 

of  Ij  will  be  lowered  by  1  (but  keep  p^ V  >  1). 
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In  general,  lowering  the  degree  could  be  dangerous.  Consider  the 
following  situation:  Let  e^  e9  be  two  local  error  indicators,  and 


is  given  a  priori.  Suppose  that  the  two  intervals  are  of  p-type 


and 


,  (v) 

e,  <  pe 
1  max 


„  .  e(v) 

e2  emax* 


By  our  algorithm,  the  degree  of  first  interval  will  be  lowered  by  1  and 
its  local  error  indicator  will  change  to  ej  >  e^.  Meanwhile  the  degree  of 
second  interval  will  be  raised  by  1  and  its  local  error  indicator  will 
change  to  e^  <  However,  it  may  happen  that  in  the  next  step  we  have 

-  _  (v+1) 

e.  -  e 
1  max 


pe 


(v+1 ) 
max 


Then  everything  will  be  back  to  the  original  state  and  the  program  will  run 
into  a  viscious  circle. 

Observe  that  the  conditions  for  this  to  happen  are 


e2  <  pe2, 


e2  < 


Pe, 


thus 


el  ,  e2 

-5-  <  P  — — 

61  S1 


e2  e2 

P  —  —  < 

e2  6I 


2  1 


where  =  e2/e2.  Therefore  p“  >  R^R2»  If  we  choose  p  as  in  (5.2), 


the  vicious  circle  can  be  prevented. 

2)  To  find  a  good  local  error  indicator  is  not  a  simple  problem, 


since  the  h-p  version  program  produces  large  elements  and  high  degrees 


We  are  using  the  following  method: 

On  the  interval  I  (which  has  degree  p) ,  the  finite  element 
solution  can  be  written  as 


u  (x) 
s 


COql!o(x)  +  ClqM(x)  +  C2q2I)(x)  +  *' 


+  Cpq<n(x) 


where  q^g(x)  and  q|^^(x)  are  linear  such  that 


q^0(xi-l) 


-  1, 


qL^0(xi) 


=  0 


q(M(xi-i)  =  °* 


q(i^(xi)  = 


and  q^(x)  (i  >  2)  are  integrals  of  Legendre  polynomials  (transformed 
to  I  =  [x^pX.D,  We  predict  Cp+^  by  solving  a  local  stiffness  matrix, 
and  the  local  error  indicators  are  defined  to  be 


e(I,p)  = 


<T  '"p-i’™1 


lie  q 
P  P 


<i>.2m1/2 


“e)} 


There  is  no  theoretical  analysis  available  which  shows  how  good  is  this 
error  indicator.  Our  numerical  computation  shows  that  for  our  examples  it 
performes  well,  but  on  the  interval  with  singularity  in  it,  this  error 
indicator  is  low  quality. 


6.  THE  PERFORMANCE  OF  VARIOUS  VERSIONS  OF  THE  FEM,  THE  CONCLUSIONS 

In  this  section  we  will  compare  and  summarize  the  performances  of  all 
versions  of  the  finite  element  method  in  a  concrete  setting  of  an  example. 
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Let  us  consider  the  problem 


-u"  =  f 

u(0)  =  u(l)  =  0 

with  the  exact  solution 

u(x)  =  xa  -  x,  a  =  1.7 

having  a  relatively  strong  singularity  at  the  origin. 

As  before,  we  are  interested  in  the  performance  measured  by  the 
energy  norm  of  the  error.  The  graphs  plotted  in  the  double  logarithmic 
scale  will  show  the  dependence  of  the  error  on  the  number  N  of  degrees  of 
f  reedom. 

Fig.  6.1  shows  the  performance  of  the  h-version  (p  =  1,2)  for  the 
uniform,  the  optimal  radical  and  feedback,  h-version.  For  comparison,  we 
also  show  the  performance  of  the  optimal  h-p  version  (i.e.,  the  geometric 


10  100  1000 
N 


Figure  6.1. 
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_  2 

mesh  with  the  ratio  g  =  (/2  -  1)“  *  .1715  combined  with  a  linear  slope  of 
degrees  s  =  2ot  -  1  =  .4). 

The  figure  clearly  indicates  that  the  h-version  with  uniform  mesh  is 
not  acceptable.  The  h-version  with  optimal  mesh  performs  relatively  well 
but  strong  refinement  could  cause  round-off  problems.  For  p  =  2  the 
relative  accuracy  of  1%  is  achieved  with  N  a  40  and  the  ratio  of  the 
sizes  of  the  maximal  and  minimal  elements  is  10^.  The  h-p  version 
requires  N  a  35,  maximal  degree  5  and  the  ratio  of  the  series  of 

9 

elements  10  for  achieving  the  same  accuracy  of  IX.  Fig.  6.1  also  shows 

the  performance  of  the  feedback  h-version  for  p  =  1  and  p  =  2.  The 

feedback  method  is  here  adaptive  with  respect  to  the  rate  of  convergence 
u.  The  rate  is  the  same  as  the  rate  of  the  h-version  with  optimal 
mesh.  (The  feedback  approach  is  more  expensive  than  the  computation  with 
a-priori  given  radical  mesh.  Nevertheless  the  cost  is  not  too  high.)  The 
figure  shows  clearly  that  when  higher  accuracy  is  required,  then  the  dif¬ 
ference  between  the  performances  of  various  versions  increases. 

Figure  6.2  compares  the  performance  of  the  h,  p  and  h-p  versions. 
It  shows  the  performance  of  the  p-version  (uniform  p  <  10  with  the  geo¬ 
metric  mesh  (q  =  0.15)  and  m  =  2,5,10  elements.  We  also  show  the  per¬ 
formance  of  the  p-version  (with  the  same  number  of  elements  m)  when  the 
degrees  p  are  chosen  in  a  feedback  way.  For  m  =  1  and  m  =  5  the 

p-version  fails  to  achieve  accuracy  of  1%.  If  m  =  10  then  the  size  of  the 

—8 

smallest  element  is  of  order  10  and  the  accuracy  of  1%  is  achieved 
for  p  =  4.  We  see  that  here  (i.e.,  for  m  =  10)  the  p-version  performs 
in  the  certain  range  of  accuracy  similarly  as  the  hp-version.  This 
clearly  indicates  the  importance  of  the  selection  of  a  proper  mesh.  We 
also  see  here  the  typical  shape  of  the  curve  (the  S-shape)  when  in  the 


higher  accuracy  is  required. 
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In  summary  we  conclude 

1)  A  uniform  mesh  cannot  produce  accurate  results  for  a  reasonable 
cost  if  the  solution  has  a  singular  behavior. 

2)  The  proper  selection  of  the  mesh  is  essential  for  the  perform¬ 
ance.  The  mesh  can  be  constructed  a-priori  if  the  structure  of  the  solu¬ 
tion  is  known  or  it  can  be  constructed  in  a  feedback  way.  The  under¬ 
refinement  of  the  mesh  has  to  be  avoided.  The  overref ineraent  does  not 

influence  too  negatively  the  performance. 

3)  The  higher  degree  elements  with  properly  designed  mesh  perform 
better  than  elements  of  low  degrees  for  both  smooth  and  nonsmooth  solu¬ 
tion.  If  high  accuracy  is  required  then,  especially,  the  high  order 


element  perform  well. 
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4)  The  p-version  is  in  practical  computations  effective  provided 
the  mesh  is  properly  designed  and  the  required  accuracy  is  achieved  in  the 
exponential  phase.  If  the  mesh  is  not  properly  designed  then  the 
p-version  does  not  perform  well  for  singular  solution  although  better  than 
the  h-version  with  uniform  mesh. 

5)  The  feedback  methods  can  be  designed  so  that  they  are  adaptive 
with  respect  to  the  convergence  and  to  the  convergence  rate  u.  They  per¬ 
form  comparably  as  the  optimal  meshes. 

Although  our  conclusions  are  based  on  the  one  dimensional  case,  our 
results  and  computational  experience  related  to  the  two  dimensional  pro¬ 
blems  indicate  that  the  conclusion  are  valid  also  in  two  dimension  case. 

Let  us  mention  that  we  did  not  address  various  aspects  of  computa¬ 
tional  complexity  as  number  of  operations,  data  flow  problems,  etc.  These 
aspects  will  be  addressed  in  detail  in  [12]  in  the  two  dimensional  setting. 
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